# Importando as bibliotecas
import numpy as np
import matplotlib.pyplot as plt
import math
# Variáveis iniciais
e0 = 8.85e-12 # Permissividade no vácuo
Vi = 0 # Potencial da placa em z = 0
V0 = 10 # Potencial da placa em z = d
d = 0.002 # A distancia entre as placas condutoras = 2 mm
L = 0.1 # O Lado = 10 cm
N = 20 # Número de segmentos
delta = L/N # Tamanho do lado de cada elemento
# Cada lado é divido em N segmentos, totalizando N² elementos
# Porém, não iremos consiserar 1/4 faltante da placa, com isso teremos 3*((N²)/4) de cada placa
qtdElementos = int(3*((N**2)/2)) #quantidade de discretizações das duas placas serão 3*((N²)/4) + 3*((N²)/4) = 3*((N²)/2)
#Estes valores das constantes serão para a questão 2, porém iremos variar o N.
# Matriz para as coordenadas do centro de cada elemento (x, y, z)
# A matriz vai ter 2 * N² elementos (duas placas condutoras)
def coordenadasCentro(n):
matrizPosicao = []
for k in range(2):
for i in range(n):
for j in range(n):
# checamos se no eixo y <= L/2 para discretizarmos todas as células delta,
# caso não, só discretizamos com y > L/2, apenas com x <= L/2
val = [(j * delta) + (delta/2), (i * delta) + (delta/2), k * d]
if i >= (n/2) and j >= (n/2):
continue
else:
matrizPosicao.append(val)
return matrizPosicao
matrizPosicao20 = coordenadasCentro(N)
#Printando matriz
# for vet in matrizPosicao20:
# for val in vet:
# print(round(val, 3), end = ' ')
# print()
# Matriz de impedância
def impendancia(m, n, R):
# R = matrizPosicao
if m != n:
# 1/(4 * pi * e0) * delta²/sqrt[(xm - xn)² + (ym - yn)² + (zm - zn)²]
return 1/(4 * math.pi * e0) * (delta**2)/math.sqrt(((R[m][0] - R[n][0])**2) + ((R[m][1] - R[n][1])**2) + ((R[m][2] - R[n][2])**2))
# [delta/(pi * e0)] * ln(1 + sqrt(2))
return delta/(math.pi * e0) * math.log(1 + math.sqrt(2))
def matrizImpedancia(n, matrizPosicao):
Z = []
# tamanho = dimensão da placa L de baixo (Z=0) que é (3/4) * n² + dimensão da placa L de cima (z=d) que é (3/4) * n²
for i in range(n):
Z.append([])
for j in range(n):
Z[i].append(impendancia(i, j, matrizPosicao))
return Z
# Matriz de impedância para o caso inicializado (N == 20)
Z20 = matrizImpedancia(qtdElementos, matrizPosicao20)
#Printando a matriz
# for i in range(qtdElementos):
# for j in range(qtdElementos):
# print(round(Z20[i][j], 3), end = ' ')
# print()
# Matriz de tensão (potencial)
def matrizTensao(n):
V = []
# // é divisão por inteiro
for i in range(n//2): # Placa condutora em z = 0
V.append(Vi)
for i in range(n//2, n): # Placa condutora em z = d
V.append(V0)
return V
# Matriz de tensão para o caso inicializado (N == 20)
V20 = matrizTensao(qtdElementos)
# Printando a matriz
for i in range(qtdElementos):
print(V20[i], end = ' ')
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
def amplitude(imp, volt):
# Resolvendo o sistema para determinar as amplitudes dos pulsos
return np.linalg.solve(imp, volt)
# Cada elemento de A20 terá uma solução para cada parte dividida da placa condutora (que está em R)
# Amplitude dos pulsos para o caso inicializado (N == 20)
A20 = amplitude(Z20, V20)
# print(A20)
def distSuperf(amp):
# A distribuição superficial de carga pode ser aproximada por uma soma dos elementos do vetor A
distribuicaoSuperficial = sum(amp)
print("ρ =", distribuicaoSuperficial) # rho
return distribuicaoSuperficial
# Distribuição superficial de carga para o caso inicializado (N == 20)
rho20 = distSuperf(A20)
ρ = 7.387195695632595e-07
# Resultado, gráfico da amplitude dos pulsos
def ampChart(amp):
plt.title('Pulsos', fontsize=20)
plt.xlabel('Elemento', fontsize=20) # Elemento = parte da placa dividida pelo segmento (N)
plt.ylabel('Amplitude', fontsize=20) # Valor esperado
plt.plot(amp, color='blue', linewidth='2')
plt.show()
# Gráfico plotado para o caso inicializado (N == 20)
ampChart(A20)
# Gráfico 3D (distribuição de carga)
def ampChart3D(centros, amp, texto = 'Distribuição superficial de carga'):
newR = np.array(centros)
# Gráfico para z = 0
# Primeira metade dos arrays
X = newR[:len(newR)//2, 0]
Y = newR[:len(newR)//2, 1]
Z = np.array(amp[:len(amp)//2])
figr = plt.figure(figsize=plt.figaspect(0.5));
# ax = plt.axes(projection='3d');
ax = figr.add_subplot(1,2,1, projection = '3d')
c = X+Y+Z
ax.scatter(X, Y, Z, c=c, cmap='viridis', edgecolor='none');
ax.set_title('Distribuição de carga na placa com z = 0');
# plt.show();
# plt.figure();
# ax = plt.axes(projection='3d');
ax = figr.add_subplot(1,2,2, projection = '3d')
ax.plot_trisurf(X, Y, Z, cmap='viridis', edgecolor='none');
ax.set_title('Distribuição de carga na placa com z = 0');
plt.show();
print('z = d')
print()
# Gráfico para z = d
# A outra metade dos arrays
Xd = newR[len(newR)//2:len(newR), 0]
Yd = newR[len(newR)//2:len(newR), 1]
Zd = np.array(amp[len(amp)//2:len(amp)])
# plt.figure();
# ax = plt.axes(projection='3d');
figr = plt.figure(figsize=plt.figaspect(0.5));
ax = figr.add_subplot(1,2,1, projection = '3d')
c= Xd+Yd+Zd
ax.scatter(Xd, Yd, Zd, c=c, cmap='viridis', edgecolor='none');
ax.set_title('Distribuição de carga na placa com z = d');
# plt.show();
# plt.figure();
# ax = plt.axes(projection='3d');
ax = figr.add_subplot(1,2,2, projection = '3d')
ax.plot_trisurf(Xd, Yd, Zd, cmap='viridis', edgecolor='none');
ax.set_title('Distribuição de carga na placa com z = d');
plt.show();
# Gráfico 3D plotado para o caso inicializado (N == 20)
ampChart3D(matrizPosicao20, A20)
z = d
# Teste para vários N
listaN = [4, 8, 12, 20, 30, 40, 50, 60]
for k in listaN:
delta = float(L)/k # lado de cada elemento
qtdElementosK = int(3*((k**2)/2))
Rk = coordenadasCentro(k)
Zk = matrizImpedancia(qtdElementosK, Rk)
Vk = matrizTensao(qtdElementosK)
Ak = amplitude(Zk, Vk)
print('Para N =', k)
ampChart3D(Rk, Ak, k)
Para N = 4
z = d
Para N = 8
z = d
Para N = 12
z = d
Para N = 20
z = d
Para N = 30
z = d
Para N = 40
z = d
Para N = 50
z = d
Para N = 60
z = d
# Cálculo da capacitância
def capacitancia(amp, n):
# Não tem capacitância em z = 0 (tensão nula)
return sum(amp[(n//2):n]) * delta**2/V0
capacitancias = []
alcance = range(5, 80, 5)
# 5, 10, ..., 70, 75
for k in alcance:
delta = L/k # lado de cada elemento
qtdElementosK = int(3*((k**2)/2))
Rk = coordenadasCentro(k)
Zk = matrizImpedancia(qtdElementosK, Rk)
Vk = matrizTensao(qtdElementosK)
Ak = amplitude(Zk, Vk)
cap = capacitancia(Ak, qtdElementosK)
capacitancias.append(cap)
print('Capacitância para N =', k, 'é:', cap);
Capacitância para N = 5 é: -8.330586333689967e-13 Capacitância para N = 10 é: -2.989551097455612e-11 Capacitância para N = 15 é: 1.280470154347472e-10 Capacitância para N = 20 é: 5.464664783592663e-11 Capacitância para N = 25 é: 4.4924343508110533e-11 Capacitância para N = 30 é: 4.213417690321582e-11 Capacitância para N = 35 é: 4.0846181158309515e-11 Capacitância para N = 40 é: 4.0027663315589416e-11 Capacitância para N = 45 é: 3.983845572859667e-11 Capacitância para N = 50 é: 3.9292233466373106e-11 Capacitância para N = 55 é: 3.942084370155707e-11 Capacitância para N = 60 é: 3.891402877684328e-11 Capacitância para N = 65 é: 3.918459688830463e-11 Capacitância para N = 70 é: 3.867327666328806e-11 Capacitância para N = 75 é: 3.9026873227979175e-11
O objetivo é o erro (diferença entre as capacitâncias do analítico e computacional) ser próximo de 0.
plt.plot(np.array(alcance), np.array(capacitancias) - e0 * (3*(L**2)/4)/d)
plt.plot(np.array(alcance), np.zeros(len(alcance)))
plt.xlabel('N')
plt.ylabel('Capacitância')
plt.show()
# Valor analítico aproximado
print('Capacitância esperada =', e0 * (3*(L**2)/4)/d)
Capacitância esperada = 3.318750000000001e-11